if (!require(tidyverse)) { install.packages('tidyverse') }; require(tidyverse)
if (!require(here)) { install.packages('here') }; require(here)
if (!require(raster)) { install.packages('raster') }; require(raster)
if (!require(sf)) { install.packages('sf') }; require(sf)
if (!require(rasterVis)) { install.packages('rasterVis') }; require(rasterVis)

First we need to to load our raster image. This is a multiband raster so we use the stack function from the raster package. For a singleband raster we would just use the raster function.

img <- stack(here('Data', '20210804_162603_13_2419_3B_AnalyticMS_SR_clip.tif'))

Now we can check the raster for some important information.

res(img) # pixel resolution
crs(img) # coordinate reference system (and units of resolution)
names(img) # band names (default)
## [1] 3 3
## CRS arguments:
##  +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## [1] "X20210804_162603_13_2419_3B_AnalyticMS_SR_clip.1"
## [2] "X20210804_162603_13_2419_3B_AnalyticMS_SR_clip.2"
## [3] "X20210804_162603_13_2419_3B_AnalyticMS_SR_clip.3"
## [4] "X20210804_162603_13_2419_3B_AnalyticMS_SR_clip.4"

For ease of reference we can give the bands better names at this point.

names(img) <- c('R', 'G', 'B', 'NIR')
names(img)
## [1] "R"   "G"   "B"   "NIR"

Now using plotRGB we can visualize the image.

plotRGB(img, 1, 2, 3)
## Error in grDevices::rgb(RGB[, 1], RGB[, 2], RGB[, 3], alpha = alpha, maxColorValue = scale): color intensity 324, not in 0:255

Well that didn’t work. The plotRGB function assumes we are using pixel brightness values (16-bit) so it struggles with our surface reflectance values being 256-bit. Here are a few ways to fix that issue.

plotRGB(img, 1, 2, 3, stretch = 'hist') # Histogram stretch

plotRGB(img, 1, 2, 3, stretch = 'lin') # Linear (scaling) stretch

# Both of these will also need stretches applied, but its good to know the scale argument
plotRGB(img, 1, 2, 3, scale = 256^2) # Theoretical maximum value

plotRGB(img, 1, 2, 3, scale = 13108) # True maximum value

Okay, at 3x3 meters, these images are slow to render. To get an extent object for subsetting you could use the drawExtent function. Obviously I can’t do that in a non-interactive session so I’m going to build one from scratch.

aoi <- c(733787.6, 4437162.0, 741327.3, 4443901.0) %>% 
  matrix(ncol = 2) %>% 
  extent()

plotRGB(img, 1, 2, 3, stretch = 'lin', ext = aoi)

To speed things up even more we can make a smaller raster by clipping.

img2 <- img %>% raster::crop(aoi)

plotRGB(img2, 1, 2, 3, stretch='lin')

Raster math is pretty straightforward.

ndvi_calc <- . %>% {(.$NIR-.$R)/(.$NIR+.$R)}

ndvi <- img2 %>% ndvi_calc

plot(ndvi)

Now we can plot NDVI as raw values, or set up a logical statement to isolate some values. Now would be a perfect time to use levelplot.

levelplot(ndvi, cuts = 5, margin=list(draw = F))

levelplot(ndvi > 0.85, cuts=1, margin=list(draw = F), colorkey=F)

We can always plot NDVI alone, but sometimes its better to combine the stacks and simply assign it to a color band (green in this case).

img3 <- stack(img2, ndvi)
plotRGB(img3, 1, 5, 3, stretch = 'lin')

Sample and Extract Values